clear
tic
syms x y z phi R L h G E mu psi rho K M
syms psi_eta psi_omega psi_xi 
b1=1.875/L;
b2=4.694/L;
r1=(sin(b1*L)-sinh(b1*L))/(cos(b1*L)+cosh(b1*L));
f1=cos(b1*x)-cosh(b1*x)+r1*(sin(b1*x)-sinh(b1*x));
r2=(sin(b2*L)-sinh(b2*L))/(cos(b2*L)+cosh(b2*L));
f2=cos(b2*x)-cosh(b2*x)+r2*(sin(b2*x)-sinh(b2*x));

nd=2; % �ھ���
% �����λ�ƺ���
% eta_base=[[sin(2*phi),sin(3*phi),sin(4*phi)].*f1,[sin(2*phi),sin(3*phi),sin(4*phi)].*f2];   % ����λ��
% omega_base=[[cos(2*phi),cos(3*phi),cos(4*phi)].*f1,[cos(2*phi),cos(3*phi),cos(4*phi)].*f2];  %����λ�ƣ��Ӷȣ�
% xi_base=[[cos(2*phi),cos(3*phi),cos(4*phi)].*sin(pi*x/2/L),[cos(2*phi),cos(3*phi),cos(4*phi)].*sin(3*pi*x/2/L)];  %����λ��
eta_base=[sin(nd*phi).*x^2,sin(nd*phi).*x^3,sin(nd*phi).*x^4];   % ����λ��
omega_base=[cos(nd*phi).*x^2,cos(nd*phi).*x^3,cos(nd*phi).*x^4];  %����λ�ƣ��Ӷȣ�
xi_base=[cos(nd*phi).*sin(pi*x/2/L),cos(nd*phi).*sin(3*pi*x/2/L)];  %����λ��

len_eta=length(eta_base);
len_omega=length(omega_base);
len_xi=length(xi_base);

% �����������꣨����ʱ��ĺ�����
for i=1:len_eta
    psi_eta(i,1)=sym(['psi_eta',num2str(i)]);
end

for i=1:len_omega
    psi_omega(i,1)=sym(['psi_omega',num2str(i)]);
end

for i=1:len_xi
    psi_xi(i,1)=sym(['psi_xi',num2str(i)]);
end

psi=[psi_eta;psi_omega;psi_xi];
psi_len=length(psi);

% �������������µ�������λ��
eta=vpa(eta_base*psi_eta);
omega=vpa(omega_base*psi_omega);
xi=vpa(xi_base*psi_xi);

% ����������Ҫ��΢��
d_omega_d_x2=diff(omega,x,2);
d_omega_d_phi2=diff(omega,phi,2);
d_omega_d_x_d_phi=diff(diff(omega,x,1),phi,1);

D0=E*h^3/12/(1-mu^2);  %��������ն�

%��������
D1=simplify(E/(1-mu^2)*[1,mu,0;mu,1,0;0,0,(1-mu)/2]);
eps_x=diff(xi,x,1)-z*d_omega_d_x2;
eps_phi=diff(eta,phi,1)/R+omega/R-z/R^2*d_omega_d_phi2;
eps_x_phi=diff(eta,x,1)+diff(xi,phi,1)/R-2*z/R*d_omega_d_x_d_phi;
eps=[eps_x;eps_phi;eps_x_phi];
% U_tmp=int(eps.'*D1*eps,z,-0.5*h,0.5*h)*R;
% U_tmp=d_omega_d_x2^2+2*mu*d_omega_d_x2*d_omega_d_phi2/R^2+(d_omega_d_phi2/R^2)^2+2*(1-mu)*(d_omega_d_x_d_phi/R)^2*R;
% U_all=0.5*int(int(U_tmp,phi,0,2*pi),x,0,L);
U_tmp=int(eps_x^2+eps_phi^2+2*mu*eps_x*eps_phi+(1-mu)/2*eps_x_phi^2,z,-0.5*h,0.5*h);
U_all=E*R/2/(1-mu^2)*int(int(U_tmp,phi,0,2*pi),x,0,L);

% ����նȾ���
for i=1:psi_len
    for j=1:psi_len
        K(i,j)=diff(diff(U_all,psi(i),1),psi(j),1);
    end
end

% ���㶯��
T_tmp=eta^2+omega^2+xi^2;
T_all=0.5*rho*h*R*int(int(T_tmp,phi,0,2*pi),x,0,L);


% ������������
for i=1:psi_len
    for j=1:psi_len
        M(i,j)=diff(diff(T_all,psi(i),1),psi(j),1);
    end
end
K_fin=double(subs(K,[R L h G E mu rho],[50 100 1 76923 200000 0.3 8e-9]));
M_fin=double(subs(M,[R L h G E mu rho],[50 100 1 76923 200000 0.3 8e-9]));
[V,D]=eig(K_fin,M_fin);
disp(V);
disp(sqrt(D));
toc

% ����ĳһ���͵�ͼ��
R_=50;
L_=100;
order=1;  % ѡ��һ������
rate=0.03;
fun_eta=matlabFunction(subs(eta,[psi;L],[V(:,order);L_]),'Vars',[x,phi]);
fun_xi=matlabFunction(subs(xi,[psi;L],[V(:,order);L_]),'Vars',[x,phi]);
fun_omega=matlabFunction(subs(omega,[psi;L],[V(:,order);L_]),'Vars',[x,phi]);
fun_x=@(x,phi)x+fun_xi(x,phi)*rate;
fun_y=@(x,phi)(R_+fun_omega(x,phi)*rate)*cos(phi)-fun_xi(x,phi)*sin(phi)*rate;
fun_z=@(x,phi)(R_+fun_omega(x,phi)*rate)*sin(phi)+fun_xi(x,phi)*cos(phi)*rate;

fsurf(fun_x,fun_y,fun_z,[0,L_,0,2*pi])
axis equal